In this tutorial we will see how to calculate the $t_0$ with the known sample. To ensure accurate estimation you should select your peak within the big mc range. You should also make your calculation for a small volume specially and for a small variation of voltage. Another pont is that you should choose a narrow window around your peaks.
# Activate intractive functionality of matplotlib
%matplotlib ipympl
# Activate auto reload
%load_ext autoreload
%autoreload 2
%reload_ext autoreload
# import libraries
import os
import numpy as np
from sklearn import linear_model
import matplotlib.pyplot as plt
from ipywidgets import fixed, interact_manual, widgets
from ipywidgets import Output
from IPython.display import clear_output
# Local module and scripts
from pyccapt.calibration.calibration_tools import share_variables
from pyccapt.calibration.data_tools import data_tools, data_loadcrop, dataset_path_qt
from pyccapt.calibration.mc import mc_tools
from pyccapt.calibration.calibration_tools import widgets as wd
from pyccapt.calibration.calibration_tools import mc_plot
By clicking on the button below, you can select the dataset file you want to use. The dataset file can be in various formats, including HDF5, EPOS, POS, ATO, and CSV.
button = widgets.Button(
description='load dataset',
)
@button.on_click
def open_file_on_click(b):
"""
Event handler for button click event.
Prompts the user to select a dataset file and stores the selected file path in the global variable dataset_path.
"""
global dataset_path
dataset_path = dataset_path_qt.gui_fname().decode('ASCII')
button
!conda install --yes --prefix {sys.prefix} pytables
tdc, pulse_mode, flightPathLength_d, t0_d, max_mc, det_diam = wd.dataset_instrument_specification_selection()
display(tdc, flightPathLength_d, max_mc)
# exctract needed data from Pandas data frame as an numpy array
# create an instance of the Variables opject
variables = share_variables.Variables()
variables.pulse_mode = pulse_mode
dataset_main_path = os.path.dirname(dataset_path)
dataset_main_path = os.path.dirname(dataset_main_path)
dataset_name_with_extention = os.path.basename(dataset_path)
variables.dataset_name = os.path.splitext(dataset_name_with_extention)[0]
variables.result_data_path = dataset_main_path + '/t_0_flight_path_distance/'
variables.result_data_name = variables.dataset_name
variables.result_path = dataset_main_path + '/t_0_flight_path_distance/'
if not os.path.isdir(variables.result_path):
os.makedirs(variables.result_path, mode=0o777, exist_ok=True)
# Create data farame out of hdf5 file dataset
data = data_tools.load_data(dataset_path, tdc.value, mode='processed')
# extract data from the path and create the Variable object
data_tools.extract_data(data, variables, flightPathLength_d.value, max_mc.value)
The maximum time of flight: 5010
data
| x (nm) | y (nm) | z (nm) | mc_c (Da) | mc (Da) | high_voltage (V) | pulse | start_counter | t_c (ns) | t (ns) | x_det (cm) | y_det (cm) | pulse_pi | ion_pp | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 0.0 | 0.0 | 0.0 | 0.0 | 14.136714 | 5019.720215 | 1003.943970 | 3495 | 0.0 | 446.853577 | 2.964898 | -0.169796 | 0 | 0 |
| 1 | 0.0 | 0.0 | 0.0 | 0.0 | 29.616535 | 5019.720215 | 1003.943970 | 3565 | 0.0 | 616.451904 | -1.936327 | 0.088163 | 70 | 2 |
| 2 | 0.0 | 0.0 | 0.0 | 0.0 | 30.456534 | 5019.720215 | 1003.943970 | 4103 | 0.0 | 623.172729 | -1.648980 | 0.672653 | 538 | 1 |
| 3 | 0.0 | 0.0 | 0.0 | 0.0 | 29.550233 | 5019.720215 | 1003.943970 | 4134 | 0.0 | 616.253052 | -1.975510 | -0.218776 | 31 | 1 |
| 4 | 0.0 | 0.0 | 0.0 | 0.0 | 28.966265 | 5019.720215 | 1003.943970 | 4205 | 0.0 | 605.431091 | 1.296327 | -0.173061 | 71 | 1 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 10201104 | 0.0 | 0.0 | 0.0 | 0.0 | 141.314648 | 6347.270020 | 1269.453979 | 8768 | 0.0 | 1154.633423 | -0.734694 | -1.577143 | 403 | 1 |
| 10201105 | 0.0 | 0.0 | 0.0 | 0.0 | 29.461296 | 6347.270020 | 1269.453979 | 8799 | 0.0 | 548.825195 | 0.408163 | 1.508571 | 31 | 1 |
| 10201106 | 0.0 | 0.0 | 0.0 | 0.0 | 29.600126 | 6347.270020 | 1269.453979 | 9443 | 0.0 | 547.803345 | -1.165714 | 0.097959 | 485 | 1 |
| 10201107 | 0.0 | 0.0 | 0.0 | 0.0 | 29.719453 | 6347.270020 | 1269.453979 | 9771 | 0.0 | 555.038513 | -1.645714 | 1.296327 | 328 | 1 |
| 10201108 | 0.0 | 0.0 | 0.0 | 0.0 | 29.823040 | 6347.270020 | 1269.453979 | 10327 | 0.0 | 556.574707 | -0.982857 | 1.933061 | 556 | 1 |
10201109 rows × 14 columns
Extract the data from the dataset file in numpy array format.
# exctract needed data from Pandas data frame as an numpy array
# exctract needed data from Pandas data frame as an numpy array
pulse_mode.value = 'voltage'
variables.dld_high_voltage = data['high_voltage (V)'].to_numpy()
variables.dld_pulse = data['pulse'].to_numpy()
variables.dld_t = data['t (ns)'].to_numpy()
variables.dld_x_det = data['x_det (cm)'].to_numpy()
variables.dld_y_det = data['y_det (cm)'].to_numpy()
# calculate the mc based on t_0 = 0
variables.mc = mc_tools.tof2mc(variables.dld_t, 0, variables.dld_high_voltage, variables.dld_x_det, variables.dld_y_det, flightPathLength_d.value,
variables.dld_pulse, mode=pulse_mode.value)
interact_manual(data_loadcrop.plot_crop_experiment_history, data=fixed(data), variables=fixed(variables), max_tof=widgets.FloatText(value=variables.max_tof), frac=widgets.FloatText(value=1.0),
bins=fixed((1200,800)), figure_size=fixed((7,3)),
draw_rect=fixed(False), data_crop=fixed(False), pulse_plot=widgets.Dropdown(options=[('False', False), ('True', True)]), dc_plot=widgets.Dropdown(options=[('True', True), ('False', False)]),
pulse_mode=widgets.Dropdown(options=[('voltage', 'voltage'), ('laser', 'laser')]), save=widgets.Dropdown(options=[('True', True), ('False', False)]),
figname=widgets.Text(value='exp_hist'));
C:\Users\APTUser\AppData\Local\Temp\ipykernel_10944\2178511811.py:1: DeprecationWarning: on_submit is deprecated. Instead, set the .continuous_update attribute to False and observe the value changing with: mywidget.observe(callback, 'value'). interact_manual(data_loadcrop.plot_crop_experiment_history, data=fixed(data), variables=fixed(variables), max_tof=widgets.FloatText(value=variables.max_tof), frac=widgets.FloatText(value=1.0),
plot the histogram of the mass spectrum.
interact_manual(
mc_plot.hist_plot,
variables=fixed(variables),
bin_size=widgets.FloatText(value=0.1),
log=widgets.Dropdown(options=[('True', True), ('False', False)]),
target=widgets.Dropdown(options=[('mc', 'mc'), ('tof', 'tof')]),
mode=widgets.Dropdown(options=[('normal', 'normal'), ('normalized', 'normalized')]),
prominence=widgets.IntText(value=100),
distance=widgets.IntText(value=100),
percent=widgets.IntText(value=50),
selector=fixed('peak'),
figname=widgets.Text(value='hist'),
lim=widgets.IntText(value=variables.max_mc),
peaks_find_plot=widgets.Dropdown(options=[('True', True), ('False', False)]),
peaks_find=fixed(True),
range_plot=fixed(False),
plot_ranged_ions=fixed(False),
ranging_mode=fixed(False),
selected_area_specially=fixed(False),
selected_area_temporally=fixed(False),
save_fig=widgets.Dropdown(options=[('True', True), ('False', False)]),
print_info=fixed(True),
figure_size=fixed((9, 5)));
C:\Users\APTUser\AppData\Local\Temp\ipykernel_10944\966379802.py:1: DeprecationWarning: on_submit is deprecated. Instead, set the .continuous_update attribute to False and observe the value changing with: mywidget.observe(callback, 'value'). interact_manual(
For the selected peak you have to select the mass-to-charge ratio of the peak. You can select the mass-to-charge ratio from the dropdown list below. You can also select the charge of the peak. After that you can add the peak to the list by clicking on the add button. You can also delete the selected peak from the list by clicking on the delete button. You can reset the list by clicking on the reset button.
isotopeTableFile = '../../../files/isotopeTable.h5'
dataframe = data_tools.read_hdf5_through_pandas(isotopeTableFile)
elementsList = dataframe['element']
elementIsotopeList = dataframe['isotope']
elementMassList = dataframe['weight']
abundanceList = dataframe['abundance']
elements = list(zip(elementsList, elementIsotopeList, elementMassList, abundanceList))
dropdownList = []
for element in elements:
tupleElement = ("{} ({}) ({:.2f})".format(element[0], element[1], element[3]), "{}({})[{}]".format(element[0], element[1], element[2]))
dropdownList.append(tupleElement)
chargeList = [(1,1,),(2,2,),(3,3,),(4,4,)]
dropdown = wd.dropdownWidget(dropdownList,"Elements")
dropdown.observe(wd.on_change)
chargeDropdown = wd.dropdownWidget(chargeList,"Charge")
chargeDropdown.observe(wd.on_change_charge)
wd.compute_element_isotope_values_according_to_selected_charge()
buttonAdd = wd.buttonWidget("ADD")
buttonDelete = wd.buttonWidget("DELETE")
buttonReset = wd.buttonWidget("RESET")
display(dropdown, chargeDropdown, buttonAdd, buttonDelete, buttonReset)
def buttonAdd_f(b,):
with out:
clear_output(True)
wd.onClickAdd(b, variables)
display()
def buttonDelete_f(b,):
with out:
clear_output(True)
wd.onClickDelete(b, variables)
display()
def buttonResett_f(b,):
with out:
clear_output(True)
wd.onClickReset(b, variables)
display()
listMaterial = buttonAdd.on_click(buttonAdd_f)
buttonDelete.on_click(buttonDelete_f)
buttonReset.on_click(buttonResett_f)
# listMaterial = buttonAdd.on_click(wd.onClickAdd)
# buttonDelete.on_click(wd.onClickDelete)
# buttonReset.on_click(wd.onClickReset)
out = Output()
display(out)
# The correct peak location based on your sample
peaks_ideal = variables.list_material
print('peak ideal mass:', peaks_ideal)
peak_x = variables.peaks_x_selected
print('peak actual mass:', peak_x)
print('The peak index are:', variables.peaks_index_list)
peak ideal mass: [1.01, 13.49, 26.98] peak actual mass: [1.0024705770121665, 14.435576305489002, 29.272140841418643] The peak index are: [0, 1, 2]
Here for each peak a mask will be created.
pick_ions_plot = []
mask = np.zeros(len(variables.mc), dtype=bool)
mc_seb_ideal = np.zeros(len(variables.mc))
# creat mask for each peak base on the peak loc. and window size
index = 0
for i in variables.peaks_index_list:
mask_tmp = np.logical_and((variables.x_hist[round(variables.peak_widths[2][i])] < variables.mc), (variables.mc < variables.x_hist[round(variables.peak_widths[3][i])]))
# Count the number of True values in the original mask_tmp
true_indices = np.where(mask_tmp)[0]
num_true_values = len(true_indices)
# If there are more than 5000 True values, randomly choose 5000 of them
if num_true_values > 5000:
selected_indices = np.random.choice(true_indices, size=5000, replace=False)
# Create a new mask with the same length as the original, initialized with False
new_mask = np.full_like(mask_tmp, False)
# Set the selected indices to True in the new mask
new_mask[selected_indices] = True
# Update the original mask_tmp with the new mask
mask_tmp = new_mask
print('peak left and right sides are:', variables.x_hist[round(variables.peak_widths[2][i])], variables.x_hist[round(variables.peak_widths[3][i])])
mask = np.logical_or(mask, mask_tmp)
mc_seb_ideal[mask_tmp==True] = peaks_ideal[index]
index += 1
bb = np.zeros(len(variables.mc))
pick_ions_plot.append(mask_tmp)
peak left and right sides are: 1.0024705770121665 1.1027176346873666 peak left and right sides are: 14.134835132463401 14.836564536189803 peak left and right sides are: 28.570411437692243 29.873623187469846
Here we plot ions in each peak base on the TOF and (x,y) position
for i in range(len(pick_ions_plot)):
fig1, ax1 = plt.subplots(figsize=(6, 6))
dld_x_masked = variables.dld_x_det[pick_ions_plot[i]]
dld_y_masked = variables.dld_y_det[pick_ions_plot[i]]
dld_t_masked = variables.dld_t[pick_ions_plot[i]]
if len(dld_x_masked) > 1000:
mask_plot = np.random.randint(0, len(dld_x_masked), 1000)
x = plt.scatter(dld_x_masked[mask_plot], dld_t_masked[mask_plot], color="blue", label='X', alpha=0.1)
y = plt.scatter(dld_y_masked[mask_plot], dld_t_masked[mask_plot], color="red", label='Y', alpha=0.1)
else:
x = plt.scatter(dld_x_masked, dld_t_masked, color="blue", label='X', alpha=0.1)
y = plt.scatter(dld_y_masked, dld_t_masked, color="red", label='Y', alpha=0.1)
ax1.set_ylabel("Time of flight (ns)", color="red", fontsize=10)
ax1.set_xlabel("position (cm)", color="red", fontsize=10)
plt.grid(color='aqua', alpha=0.3, linestyle='-.', linewidth=2)
plt.legend(handles=[x, y], loc='upper right')
plt.ylim(bottom=plt.yticks()[0][0], top=plt.yticks()[0][-1])
plt.xlim(left=plt.xticks()[0][0], right=plt.xticks()[0][-1])
plt.savefig(variables.result_path + 'position_peak.png', format="png", dpi=600)
plt.savefig(variables.result_path + 'position_peak.svg', format="svg", dpi=600)
plt.show()
As you saw the TOF changes base on the (x,y) of the events. Therefore we creat a mask to only select the ions in center (8m m*8mm) of detector. This helps to cansel out the variation in TOF due to hit position.
Reformulate the equation:
# use mask_equal to have equal number of ions for each peak
# only peak the value in the center of detector 10mm * 10mm
detector_squre = 0.5
vol_variation = 5500 # peak only for 200 voltage variation for vol_variation +/- 100
fig1, ax1 = plt.subplots(figsize=(6, 6))
dld_x_masked = variables.dld_x_det[mask]
dld_y_masked = variables.dld_y_det[mask]
dld_t_masked = variables.dld_t[mask]
dld_highVoltage_masked = variables.dld_high_voltage[mask]
dld_pulse_masked = variables.dld_pulse[mask]
mask_tmp_middle = np.logical_and((np.abs(dld_x_masked) < detector_squre), (np.abs(dld_y_masked) < detector_squre))
mask_tmp_mvoltage = mask_tmp_mvoltage = np.logical_and((dld_highVoltage_masked < (np.mean(dld_highVoltage_masked)+200)), (dld_highVoltage_masked > (np.mean(dld_highVoltage_masked)-200)))
mask_tmp_mvoltage = np.ones(len(mask_tmp_middle), dtype=bool)
mask_f = np.logical_and(mask_tmp_mvoltage, mask_tmp_middle)
dld_x_masked = dld_x_masked[mask_f]
dld_y_masked = dld_y_masked[mask_f]
dld_t_masked = dld_t_masked[mask_f]
dld_highVoltage_masked = dld_highVoltage_masked[mask_f]
dld_pulse_masked = dld_pulse_masked[mask_f]
mc_seb_reg_masked = mc_seb_ideal[mask]
mc_seb_reg_masked = mc_seb_reg_masked[mask_f]
x = plt.scatter(dld_x_masked, dld_t_masked, color="blue", label='X', alpha=0.1)
y = plt.scatter(dld_y_masked, dld_t_masked, color="red", label='Y', alpha=0.1)
ax1.set_ylabel("Time of flight (ns)", color="red", fontsize=10)
ax1.set_xlabel("position (cm)", color="red", fontsize=10)
plt.grid(color='aqua', alpha=0.3, linestyle='-.', linewidth=2)
plt.legend(handles=[x, y], loc='upper right')
plt.ylim(bottom=plt.yticks()[0][0], top=plt.yticks()[0][-1])
plt.xlim(left=plt.xticks()[0][0], right=plt.xticks()[0][-1])
plt.savefig(variables.result_path + 'center.png' , format="png", dpi=600)
plt.savefig(variables.result_path + 'center.svg' , format="svg", dpi=600)
plt.show()
We calculate the ${t_0}$ base on:
seb_t = dld_t_masked * 1E-9 # tof in s
seb_factor = np.sqrt(mc_seb_reg_masked * 1.66E-27 / (2 * 1.6E-19 * dld_highVoltage_masked))
seb_factor = seb_factor * 1E6
seb_t = seb_t * 1E9
t0_seb_fixed = np.mean(np.array([seb_t]).squeeze(0) - (flightPathLength_d.value * np.array([seb_factor]).squeeze(0).reshape(-1, 1)))
print('Linear fixed path lenght -- the flight path lenght(slop): {:.2f}'.format(flightPathLength_d.value), '(mm)', '\nthe corrected t_0(intercept): {:.2f}'.format(t0_seb_fixed), '(ns)')
Linear fixed path lenght -- the flight path lenght(slop): 110.00 (mm) the corrected t_0(intercept): -13.83 (ns)
# Plot outputs
# fig1, ax1 = plt.subplots(figsize=(5.5/2.54, 5.5/2.54))
fig1, ax1 = plt.subplots(figsize=(5.5, 5.5))
peaks_data = plt.scatter(seb_factor, seb_t, label='Ions', alpha=1, color='slategray')
axes = plt.gca()
x_vals = np.array(axes.get_xlim())
linear_fix, = plt.plot(x_vals, t0_seb_fixed + flightPathLength_d.value * x_vals, '--', color='red', label='line' )
plt.grid(color='aqua', alpha=0.3, linestyle='-.', linewidth=2)
# plt.legend(handles=[linear_fix, peaks_data], loc='lower right')
ax1.set_ylabel("Time of Flight (ns)", fontsize=8)
ax1.set_xlabel(r"$\sqrt{\frac{\frac{m}{n}}{2e \alpha V}} (ns/mm)$", fontsize=8)
plt.tight_layout()
plt.ylim(bottom=plt.yticks()[0][0], top=plt.yticks()[0][-1])
plt.xlim(left=plt.xticks()[0][0], right=plt.xticks()[0][-1])
plt.savefig(variables.result_path + 'fixed_flight_path.svg', format="svg", dpi=600)
plt.savefig(variables.result_path + 'fixed_flight_path.png', format="png", dpi=600)
plt.show()
linear = linear_model.LinearRegression()
linear.fit(np.array([seb_factor]).squeeze(0).reshape(-1, 1), np.array([seb_t]).squeeze(0))
d_seb = linear.coef_.item()
t0_seb = linear.intercept_.item()
print('Linear -- the corrected flight path lenght(slop): {:.2f}'.format(d_seb), '(mm)', '\nthe corrected t_0(intercept): {:.2f}'.format(t0_seb), '(ns)')
Linear -- the corrected flight path lenght(slop): 106.98 (mm) the corrected t_0(intercept): -2.58 (ns)
# Plot outputs
# fig1, ax1 = plt.subplots(figsize=(5.5/2.54, 5.5/2.54))
fig1, ax1 = plt.subplots(figsize=(5.5, 5.5))
peaks_data = plt.scatter(seb_factor, seb_t, color="slategray", label='Ion', alpha=0.1)
axes = plt.gca()
x_vals = np.array(axes.get_xlim())
linear, = plt.plot(x_vals, t0_seb + d_seb * x_vals, '--', color='r', label='fit' )
plt.grid(color='aqua', alpha=0.3, linestyle='-.', linewidth=2)
# plt.legend(handles=[peaks_data, linear, rigid, huber,lasso], loc='lower right')
ax1.set_ylabel("Time of flight (ns)", fontsize=8)
ax1.set_xlabel(r"$\sqrt{\frac{\frac{m}{n}}{2e \alpha V}} (ns/mm)$", fontsize=8)
plt.legend(handles=[linear, peaks_data], loc='lower right')
plt.tight_layout()
plt.ylim(bottom=plt.yticks()[0][0], top=plt.yticks()[0][-1])
plt.xlim(left=plt.xticks()[0][0], right=plt.xticks()[0][-1])
plt.savefig(variables.result_path + 'regression.svg', format="svg", dpi=600)
plt.savefig(variables.result_path + 'regression.png', format="png", dpi=600)
plt.show()
Plot the m/c with new ${t_0}$:
# recalulate the mc based on the new t_0
mc = mc_tools.tof2mc(variables.dld_t, 14, variables.dld_high_voltage, variables.dld_x_det, variables.dld_y_det, flightPathLength_d.value,
variables.dld_pulse, mode=pulse_mode.value)
mc_hist = mc_plot.AptHistPlotter(mc[mc < 40], variables)
y, x = mc_hist.plot_histogram(bin_width=0.1, label='mc', steps='stepfilled', log=True, grid=True, fig_size=(9, 5))
peaks, properties, peak_widths, prominences = mc_hist.find_peaks_and_widths(prominence=100, distance=10, percent=50)
mc_hist.plot_peaks()
mc_hist.plot_hist_info_legend(label='mc', bin=0.1, background=None, loc='right')